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The Bethe-Salpeter equation and leptonic decay constants for pseudoscalar and vector quark- 
antiquark mesons with arbitrary quark masses are studied in Coulomb gauge, under a leading order 
truncation. As input, we use a pure linear rising potential, supplemented by a contact term arising 
from the conservation of total color charge. It is shown how the equations can be written in terms 
of manifestly finite functions, despite the infrared singular interaction. The resulting equations are 
solved numerically. Both the pattern of dynamical chiral symmetry breaking and the leading order 
. . - heavy quark limit are visible. 
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> \ I. INTRODUCTION 

o 

■ In the theory of strong interactions, quantum chromodynamics (QCD), there are two outstanding issues of impor- 
' tance: confinement and dynamical chiral symmetry breaking. Both of these phenomena are highly nontrivial and are 
t-h . reflected in the nonperturbative regime of the theory. Perhaps the most naive picture of confinement, at least applied 
to the quark sector of QCD, is the idea that there exists a linear rising potential between quarks at large separation. 
If a quark is pulled away from its parent hadron, the energy of the system rises until it is energetically favorable that 
a quark-antiquark pair is created and forming a new system of two hadrons whose constituents remain hidden from 
view in the same way as for the original. In momentum space, such a linear rising potential corresponds to a strongly 
infrared enhanced interaction. In Coulomb gauge, lattice results [H-[l] show that the temporal component of the 
nonperturbative gluon propagator (that mediates the quark-antiquark interaction) indeed exhibits such an infrared 
enhancement. 

The inclusion of an infrared singularity, such as that of the aforementioned temporal gluon propagator, into calcu- 
lations of finite, physical observables is not trivial. Because the singularity occurs for vanishing momentum, a natural 
place to investigate would be a nonperturbative, continuum formalism, such as that provided by the Dyson-Schwinger 
and Bethe-Salpeter equations 1 . The Coulomb gauge Dyson-Schwinger equations for QCD have been derived in both 
the first and second order formalisms (see Refs. [17H19| and references therein) and various aspects explored, such as 
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1 the one-loop perturbative behavior [18H20I ]. the Slavnov- Taylor identities [211 ] and the emergence of a nonlocal con- 
straint on the total color charge [22j . Additionally, the case of heavy quarks in the rest frame has also been studied 
[2314251 ] . These latter studies show that in the absence of pure Yang-Mills corrections, the Coulomb gauge rest frame 
£S| \ heavy quark limit (where pair creation is absent, such that the linear rising potential remains valid at very large 
separations) is described precisely by an interaction consisting of single exchange of a temporal gluon. The infrared 
singularity associated with this interaction is cancelled only when considering color singlet (meson or baryon) bound 
states, otherwise the objects (propagators or color nonsinglet states) are unphysical. 



The issue of dynamical chiral symmetry breaking was studied in Coulomb gauge many years ago (see, for example 



H . Rds. 21 i 2s] ). [1 was l'< •mid Hint wit It an infrared enhanced teniix .ml ujtton interact ion. chiral svnnnot i:v bvea.kin;-' does 
occur, although the chiral condensate and pion leptonic decay constant are too small (it has been demonstrated that 
the coupling of the quarks to transverse spatial gluons results in a significantly larger condensate [29j ) . Recently, it has 
been shown that the twin limits of the chiral and heavy quark propagator in Coulomb gauge, with an infrared enhanced 
interaction, can be accommodated within a single leading order truncation scheme for the Dyson-Schwinger equations 
[30l] 2 . The truncation scheme is centered around an Ansatz for the Coulomb kernel occurring in the action: this term 
originates from the resolution of Gauss' law and is responsible for the temporal component of the interaction between 
color charges, i.e., it is intimately connected to the temporal component of the gluon propagator. The interaction 
between quarks and the spatial components of the gluon field, the three-, and the four-gluon vertices are all neglected 
under this leading order truncation scheme. 

In this paper, the truncation scheme considered in Ref. [30f is extended to the quark-antiquark Bethe-Salpeter 
equation for pseudoscalar and vector mesons. The primary aim of this study is to show how the infrared singular 



1 An alternative approach is to consider the Coulomb gauge Hamiltonian: see for example, Refs. [9l 4l6(l and references therein. 

2 The interplay of the chiral and heavy quark symmetries has also been studied in an extended Nambu— Jona-Lasino model [3l| . 
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interaction may be consistently included into the formalism in such a way that physical quantities (meson masses and 
leptonic decay constants with arbitrary quark mass configurations) remain finite. Given that both dynamical chiral 
symmetry breaking and the heavy quark limit are qualitatively (if not quantitatively) present in the quark propagator 
under this truncation, the corresponding signals for the meson spectrum will be investigated, both analytically and 
numerically. As input, the interaction corresponding to a pure linear rising potential will be used, although it is 
recognized that this is an over-simplification from the physical standpoint. 

The paper is organized as follows. In the next section, we shall briefly review the quark propagator under the 
leading order Coulomb gauge truncation scheme. Section III then goes on to discuss the axialvector Ward-Takahashi 
identity within the context of Coulomb gauge and in the rest frame. The Bethe-Salpeter equation and leptonic decay 
constant for pseudoscalar mesons will then be discussed in Sec. IV, with particular emphasis on how the infrared 
divergence of the interaction can be rendered finite. In Sec. V, the process will be repeated for vector mesons. After 
briefly discussing the numerical implementation of the equations in Sec. VI, the results will be presented in Sec. VII. 
We finish with a summary and conclusions. 



II. QUARK GAP EQUATION 

The central input into the Bethe-Salpeter equation for quark-antiquark mesons stems from the quark propagator 
and the interaction between the quarks. Let us thus begin by briefly reviewing some pertinent results for the quark 
propagator obtained previously under the leading order Coulomb gauge truncation scheme that will form the basis 
for this study 30]. These results incorporate both dynamical chiral symmetry breaking and the heavy quark limit. 

The leading order truncation scheme is based on the following instantaneous interaction term in the Coulomb gauge 
action: 



dx dy 



(2.1) 



where p a is the color charge (superscript a denotes here the color index in the adjoint representation) associated with 
the fields. In this case we are interested only in the quark component, 

p a x = gq xl a T a q x , (2.2) 

with (conjugate) quark field at position x denoted by (q x ) q x , coupling g, and color matrix T a (Hermitian generator 
of the SU(N C ) group). In this study, we shall consider an interaction (color diagonal, and in momentum space) of the 
form 

g*C F F(q) = C(2n) 3 S(q) + (2.3) 

(q Y 

(Cf = (N% — l)/2N c is the color factor associated with the quarks). Notice that since the interaction vanishes faster 
than 1/q 2 in the ultraviolet (UV), there will be no need to consider the renormalization of any quantities and we can 
omit any discussion of the renormalization constants. The first term of the interaction arises from total color charge 
conservation (see Refs. [22|, H(| for a detailed explanation): the constant C is considered finite until the end of the 
calculation, whereupon we take the limit C — > oo. Physically, C can be thought of as an infinite shift in the potential 
connected to the fact that one requires an infinite amount of energy to create a color nonsinglet, such as a single 
quark or gluon, from the vacuum (see later). The second term is strongly infrared (IR) divergent with a coefficient 
a. The quark-antiquark potential, derived in the Coulomb gauge heavy quark limit and neglecting the contribution 
of pure Yang-Mills vertices [U [30] reads, for spatial separation r, 

V(r) = g 2 C F f dqF(q) (1 - e 19 ") = ar, (2.4) 



where dq = d<f/(27r) 3 . The coefficient a can be identified as the string tension associated with a purely linear rising 
potential. The specific value of a will not be of primary concern in this study because within the context of the 
leading order truncation scheme it will suffice that all dimensionfull quantities may be expressed in appropriate units 
of er, such that the cancellation of the infrared singularities may be demonstrated. As will be seen in Sec. VII, the 
quantitative results arising from the above potential do not lend themselves to a detailed phenomenological study 
(this would require a more sophisticated truncation scheme beyond leading order, such as in Ref. [29]). However, we 
should specify some value for a in order to make a comparison to physical quantities. The expectation value of F is 
given by the instantaneous part of the temporal gluon propagator |30l . l32j and the string tension associated with the 



3 



coefficient of the infrared divergence is called the Coulomb string tension (a c ). Thus, we should use a = a c as input. 
It is known that a c is equal to or larger than the physical Wilson string tension (<J\y) 123J and lattice results indicate 
that a c may be up to ~ 3 times lar ger than aw PH3J • At the same time though, it will also be useful to compare 
with the earlier studies of Refs. |26l - |28j where aw was used. In light of these considerations, we shall consider the 
range ^fa ~ 440 — 762MeV (the lower value corresponding to Ja~w an d the upper value to \/?>aw) to give an estimate 
for the input scale. 

Let us now turn to the quark propagator. With the interaction, Eq. (|2.1[) . and under the (leading order) truncation 
scheme to include only one loop terms, we have the following Dyson-Schwinger equation (in momentum space) for 
the quark two-point proper function, T qq [30] : 

Iw(p) = » [l°Po - TP- m] + ig 2 C F f dkF(p~ fc) 7 °W^(fc)7°, (2.5) 



where m is the bare quark mass and dk — <i 4 fc/(27r) 4 . T qq can be written in terms of two scalar dressing functions 
A p and B p (the subscript p denotes the spatial momentum dependence): 

I?,(p) = i [j°Po ~ 1-PA P - B p ] . (2.6) 
The corresponding quark propagator, W qq (p) is 

W- qq {p) = [ 7 % ~ i-pA p + B p ) , A f (p) =pl- p 2 A 2 p ~B 2 + z0+. (2.7) 

Notice that the energy dependence of both the proper two-point function and the propagator is trivial. This arises 
from the instantaneous character of the interaction and means that the dressing functions (A and B) are independent 
of the energy po. This will lead to important simplifications in the Bethe-Salpeter equation. A possible fourth Dirac 
structure, proportional to 7 °jPo7 • p and arising in the noncovariant Coulomb gauge context does not appear, just as 
in the perturbative case j 1QI ] . 

The static quark propagator is defined as 

where we define the quark mass function, M p , and quasiparticle energy, oj p : 



M p = ^, u p = y/p 2 +M*. (2.9) 

In the chiral limit, the quark condensate can be defined as (trace over Dirac matrices) 

<qq>=N c J dpTY d W^(p)^-2N c I dp^. (2.10) 

The meaning of the quark mass function becomes clear when one considers the quark Dyson-Schwinger equation in 
the context of the interaction specified in Eq. (|2.3[) . In terms of the dressing functions, A and B, Eq. (|2.5p decomposes 
to 

1 2 ^ f dkF(p- k)p-k 
7:9 C F / '— 

2 J uj k j 

1 2 f dkF(p- k) 



A p = l + -g*C F / £ '- L —, (2.11) 

£ J u k p 



B p = m + -g 2 C F / £ '-M k , (2.12) 

2 J uJk 



whereas combining the two equations in terms of the mass function leads to the Adler-Davis gap equation [27j : 

1 2ri f dkF(p-k) 



M p = m + ^9 2 C F J ■ 



p-k 
P' 



M k - "—M v 



(2.13) 



Given that the interaction F contains both the divergent coefficient C multiplying the <5-function and the IR singular 
component, A and B are both divergent quantities. However, the integrand of the gap equation contains exactly the 
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FIG. 1: [left panel] Quark mass function, M(x), and [right panel] dressing, M{x) — m, plotted as functions of x = p 2 for a 
range of quark masses. All dimensionfull quantities are in appropriate units of the string tension, a. See text for details. 



cancellation necessary to ensure that the ratio of A and B, namely, the quark mass function M, is finite. The chiral 
condensate is also finite (see also Sec. VII). Considering the quark proper two-point function written in terms of the 
self-energy, S, 

I^Op) = * [7% - TP- m q ) + £(p), (2.14) 

the divergence inherent to the dressing functions A and B has the obvious interpretation of shifting the pole of the full 
quark propagator to infinity such that one would require infinite energy to excite a single quark from the vacuum. The 
dynamical content of the quark propagator is contained within the finite mass function and we shall see that in the 
end, only the mass function (and associated quasiparticle energy) enters the Bethe-Salpeter equation for color-singlet 
mesons. 

A detailed numerical solution to Eq. (|2.13|) for a range of quark masses was presented in Ref . [30] . Two aspects will 
be of relevance to this study and for the convenience of the reader, we briefly repeat them here. The mass function 
is plotted (in appropriate units of the string tension, a) in the left panel of Fig. [JJ One can see that in the chiral 
case (to = 0) the quark mass function is nonzero, signaling the presence of dynamical chiral symmetry breaking. 
In the right panel of Fig. [TJ the dressing {M — m) is plotted. As the bare quark mass, to, increases, the dressing 
initially increases but for larger m starts to decrease. In the (heavy quark) limit to — > 00, one sees that M —> to. 
The original derivation of Eq. (|2.13[) [27] arose from considering the chiral quark limit. However, with the results for 
arbitrary quark masses it was shown [30j that the same equation reproduces the rest frame Coulomb gauge heavy 
quark limit (in the absence of pure Yang-Mills contributions) [23|. Namely, the leading order truncation considered 
here encompasses both chiral and heavy quark physics. 



III. AXIALVECTOR WARD-TAKAHASHI IDENTITY 



As is well understood, the pion can be regarded as the Goldstone boson of chiral symmetry breaking. Considering 
a chiral rotation of the quark field, 



(3.1) 



q x -> exp<; -ta )q x , 



where the a a parametrize the rotation and the r a are the Pauli matrices corresponding to a two-flavor system 
(a = 1, 2, 3), the pion vertex may be defined via the (chiral) quark self-energy as 



da a 



exp < ia b — 7 5 \ exp < ia c — 7 



{7 5 ,£(P)} 



(3.2) 



Such an approach was used to construct the pion Bethe-Salpeter equation in, for example, Ref. (26j . More generally, 
the connection between the quark gap equation and the Bethe-Salpeter equation for quark-antiquark bound states 
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is expressed via the axialvector Ward-Takahashi identity (AXWTI). In Coulomb gauge and in the chiral limit, the 
AXWTI was applied in, for example, Ref. [27|. Here, we are interested in finite, arbitrary mass quarks and the 
leptonic decay constants, so we shall closely follow the approach of Ref. [Hj]. There, the context was the consideration 
of covariant gauges. Here, for the convenience of the reader, we shall repeat their argumentation for Coulomb gauge 
in the rest frame and specifically within our leading order truncation scheme. 

Splitting the energy and spatial momentum arguments of the (noncovariant) propagators and proper functions, so 
that, for example, Tq q (p) = Tq q (p°,p), let us consider the following combination (the energy and spatial momentum 
dependence of A will be explained below): 

A Q5 (p;p°,P ) = r+ (p° + P /2,p) 7 5 ^ +7 5 yr_ ? (p° - P°/2,p), (3.3) 

where the two quark proper two-point functions, E^ , have arbitrary bare masses m ± and correspondingly, different 
dressing functions A , B and M ± . Using the quark Dyson-Schwinger equation, Eq. (|2.5[) . to rewrite the proper 
two-point functions, one obtains after some straightforward manipulation 

A a5 (p;p°,P°) = -z 7 5 [P° 7 ° + (m+ + m - )] y 



ig 2 



C F J dkF(p- fc) 7 °W+(fc° + P°/2, k)A a5 (k; k°, P°)W^ q (k° - P°/2, fc) 7 °. (3.4) 



We immediately notice that the above is a truncated Dyson-Schwinger (or inhomogeneous Bethc-Salpeter) equation 
for a generalized quark-antiquark, color-singlet vertex A a5 . Given that the right-hand side expression is independent 
of the energy p°, A a5 is also independent of p° and this arises because of the instantaneous character of the interaction. 
The energy P° is the total energy flowing through the quark-antiquark pair, whereas the spatial momentum p (or k 
within the integral) flows along the quark line. The vertex is thus in the rest frame of the quark-antiquark pair and 
this shall be extremely useful in our analysis. It is convenient to use the shorthand notation k to denote the energy 
and spatial momentum arguments (k° ± P° /2, k) (similarly for p^) such that the above equation reads 

A a5 (p;P ) = -z 7 5 [P° 7 + (m + +m-)] ll-ig 2 C F J dk F(p - k) 1 °W+ q {k + )K a5 (k- P°)Wr(fc-) 7 °. (3.5) 

The axialvector Ward-Takahashi identity tells us that this generalized vertex can be rewritten as 

A a5 {p ; P°) = -tP°Tf(p; P°) - i(m + + m~)T a5 {p; P°), (3.6) 

where Tg 5 and T a5 are the temporal component of the axialvector vertex function and the pseudoscalar vertex function, 
respectively (both in the rest frame). Separating the temporal axialvector and pseudoscalar components of A° 5 in 
Eq. (|3.5p . one can write down the Dyson-Schwinger equations for Tg 5 and r a5 : 

rg 5 (p;P°) -7 5 7°y - i9 2 C F J dkF(p- fc) 7 °W ? +(fc + )rg 5 (fc;P )M/--(fc-) 7 , (3.7) 

r a5 (p;P°) = 7 5 y - ig 2 C F J dkF(p- fc) 7 °M/+(fc+)r Q5 (fc;P°)^--(fc-) 7 . (3.8) 

Notice that both of these equations hold for arbitrary total energy P , i.e., they are not at resonance. 

Now let us consider the homogeneous Bethe-Salpeter equation in the pseudoscalar channel. Taking Eq. (I3.8P and 
making the Ansatz that T a5 has a simple pole (in the rest frame and with the quantum numbers of a pseudoscalar 
vertex) with an as yet unknown residue, rps, of the form 

rps 



ra5 (p ; p0 ) = p 2 ;, 2 r ps(p; p0 ) + no ™-> ( 3 -9) 

m PS 



(with finite parts denoted 'non-res.') then, at resonance, one has the homogeneous Bethe-Salpeter equation: 

r PS (p;P°) = PS -w 2 c F J dkF(p- fc) 7 °H^(fc + )r» s (fc;P°)w-(fc-) 7 . (3.10) 

Alternatively, taking Eq. (|3.7|) and making instead the Ansatz that 



r a o 5 (P;P°) - Z AV m2 r^(p;P") + non-res„ (3.11) 
^0 m PS 




> 
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FIG. 2: [top line] Diagrammatic representation of Eq. (|3.15p in terms of the Bethe-Salpeter kernel, K, [second line] expanded 
as a series in K and [third line] the equivalent expression in terms of the amputated four-point function G ( - 4 K [fourth line] 
The equation for G^ 4 ' (which can be also understood via the expansion in terms of K). Internal lines denote dressed quark 
propagators, the dot denotes the tree-level vertex. See text for details. 



one arrives at the identical homogeneous equation. Combining Eqs. (|3.9l3.1ip and inserting into the AXWTI, Eq. Q3.6p . 
one also has that 



A a5 (p ; P°) = -i [r AV M 2 PS + r PS (m+ + m")] ^fe^ffi + n. ,n-u<*. 



P 2 - M 2 
m PS 



(3.12) 



Given the definition of A a5 , Eq. (|3.3[) . and knowing that the quark proper two-point function contains no resonant 
components, the unknown residues rps and tav must obey the relation 



r AV Mp S 



-rps{m + + m ). 



(3.13) 



To obtain more information about the residues, some manipulation of the kernel is required. Replacing the Dirac 
(and trivial flavor) indices and denoting the Bethe-Salpeter kernel K such that 



k ot Mp, k) = ig 2 C F P(p- k) [/] [ 7 °; 



5/3 ■ 



(3.14) 



Eq. 



can be written as 



T%(p;P°) = 







7 — 
. 2 _ 


-/ 

a/3 J 



dkK ar Mp, k) w+(k+)r a5 (k-,p°)w^(k-) 



7^ 



(3.15) 



with a similar form for Eq. p.7p . Equation p.!5p is represented diagrammatically in the top line of Fig. [2j The 
equation can be recursively expanded to an infinite series in the kernel as shown in the second line of Fig. [2l One 
sees that this has the effect of replacing the nonperturbative Bethe-Salpeter vertex (r a5 ) within the integral with 
the tree-level term (j 5 T a /2) and an infinite series in terms of the kernel (K). With the form of the kernel under 
consideration here, Eq. (|3.14l) . this is nothing other than the ladder resummation of the Bethe-Salpeter equation. 
However, it is also known that the ladder resummation is directly connected to the amputated connected four-point 
quark-antiquark Green's function, (studied explicitly within the context of heavy quarks in the Coulomb gauge 
rest frame and with the same truncation scheme in Ref. [25[). The equation for is shown diagrammatically in 
the last line of Fig. [2j Replacing then the infinite ladder series with G^ in the equation for r a5 gives the third line 
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of Fig. [2j Equations (|3.7I3.8[) can thus be rewritten as 

L Z J a/3 J 



w+(k + h 5 ~/° Y w--(k-) 



#7 



a/3 



w-+(k + h 5 Y w^ q (k-) 



(3.16) 
(3.17) 
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It is explicitly known that in the Coulomb gauge heavy quark limit, G^ 4 ' contains (physical) resonant components in 
the color-singlet channel [25| . In the pseudoscalar channel under consideration here, let us write 



f*P$-M PS 



r ps (k;-P° 



7<5 



(3.18) 



where N is an unspecified normalization constant (we will use N = 1 later). The conjugate Bethe-Salpeter vertex 
function is defined with the aid of the charge conjugation matrix C, such that 



T PS (k;-P u ) = Crf s (-k;-P u )C 



(3.19) 



(T denotes the transpose). Taking the resonance Ansatz for Tg 5 , Eq. (|3 . 1 1 1) . or r a5 , Eq. (|3.9p as appropriate, along 
with Eq. (piU) for G (4 \ inserting into Eqs. q3.16j3.17j> . reorganizing the Dirac indices and evaluating the color and 
flavor traces, one obtains 



P — K'J 

r A vP° = PS -iNN c Tr d I Sk^ u W+(k + )T PS (k; P u )^_(fc-) 



Pn=Mi 



rps 



iNN c Tr d I ak^W+(k + )r PS (k;P )W^ a (k-) 



qq\ 



The identity relating the residues, Eq. (|3.13[) . can thus be written (independently of the normalization) 







Pf=Afi 



Tr d / dk^ [P ^ - (m+ +m-)]W+(k + )T PS (k;P )W q (k~). 



(3.20) 



(3.21) 



The leptonic decay constant for a pseudoscalar meson is defined as the coupling to the point axial field (see, for 
example, Refs. [26|, |3J,[3^|). In Minkowski space and in the rest frame, it can be defined as (trace over Dirac matrices) 



Nr. 



Mps 



fps = T7 £-Tr d I Sk 7 5 P° 7 W+(k+)T^ s (k;P )W^(k-) 



qq\ 



(3.22) 



where T ps is the normalized Bethe-Salpeter vertex function (the normalization will be discussed in the next section). 
It is useful to define a second quantity 



i PS =N c Tr d dk^W+(k+)T^ s (k;P^W qq (k-), 



such that Eq. (|3.2ip reads 



M Psfps = (m + + m )h PS . 
In the chiral limit, the Gell-Mann-Oakes-Renner (GMOR) relation [35[ states that 

M PS f PS m± ^° -( m + + m -) <qq> . 
Comparing with Eq. (|3. 241) . we see that in the chiral limit 



IPS 



>o <qq> 



= 2N r / dk 



PS 



PS 



(3.23) 



(3.24) 



(3.25) 



(3.26) 



showing that hps is a generalization of the chiral condensate to finite quark masses. Comparison to the GMOR 
relation will define one method of normalizing T ps aside from the canonical normalization (to be discussed later). 
Notice that with these conventions, one would have f v — 92.2MeV, a factor of \/2 smaller than the convention used 
nowadays (/„• = 130.4MeV as defined in Ref. |36j). 
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IV. PSEUDOSCALAR MESONS UNDER TRUNCATION 



Having introduced the framework for studying the Bethe-Salpeter equation in Coulomb gauge and under the leading 
order truncation, let us now show how this framework can be applied. In this section, we shall be concerned with 
the pseudoscalar mesons; in the next section we shall repeat the analysis for vector mesons. Our aim is to show 
how the homogeneous Bethe-Salpeter equation for the meson mass, Mps, the leptonic decay constant, fps, and 
the generalized condensate, hps> may be rewritten in terms of explicitly finite expressions despite the IR divergent 
interaction. 

Removing the (overall) flavor factor, the general Dirac decomposition for the normalized pseudoscalar Bethe-Salpeter 
vertex, at resonance and in the rest frame, can be written 

r$ s (p; p°) = 7 5 [r$(p 2 ) + F°7°rf (p 2 ) + rpv%{f) + PVrpif (p 2 )] . (4.i) 

In principle, the scalar dressing functions have the arguments F^(p; P°). However in practice, since the homoge- 
neous Bethe-Salpeter equation used to derive these functions is valid only at resonance (P| = M ps ), the are scalar 
functions of p alone and P° is merely a label. To condense the notation, we shall write the momentum dependence as 
a subscript, i.e., T^. Inserting the vertex decomposition, Eq. (|4.1[) . into the expressions for fps and hps, Eqs. (I3.22[) 
and (13.231) . respectively, expanding the quark propagators using Eq. (|2.7[) , performing the Dirac traces and the energy 
integrals, one arrives at the following expressions: 

dk 



fps = ^N c / r r - + ..-, 2l {Kk + M fc " W +] + r 2 ^fc 2 K+ - u- k ] 

J u k l/o - (^k +UJ k) \ L 

i^k + u k\ K 



hps = 2iN, 



[P Q 2 - (u>+ + Cj-) 2 ] \ k k [M+ + M~] V 0fc L k k k kl kl J 
-P 2 (lf fc [M+u- + M-u+] + Itlfc 2 [ u + + u-] ) } , (4.2) 



where we have introduced the notation 



(in effect, the incorporate all the information about the IR divergence of the two quark self-energies and how they 
enter the expressions). We notice that in the above, the following common combinations of functions arise: 

Qft - r& [m+ljk + M fe -w+] + r^i 2 [ u + - u>t] , 

Qik = Vik [K"k + M-ut\ + r 3 w fc fc 2 K+ + u-} , (4.4) 

allowing us to write 

[ p o - Pit +"k) 





1 ps = 2^ c / _ + __ rn2 ,_+,.■.._„■, <! (<# + fiifc ) r L a r + ; .^ QgL - P<j Qf & }> ■ (4.5) 



As emphasized, one of the most important features of the above two quantities is that both are physical (hps reduces, 
up to constant prefactors, to the chiral condensate in the chiral limit) and, in particular, cannot contain unphysical IR 
divergences such as those contained within the functions d) ± (i.e., proportional to A^). As mentioned previously, in 
the Coulomb gauge heavy quark limit it is known that the amputated connected four-point quark-antiquark Green's 
function, G^ 4 \ explicitly contains physical resonance poles [25| . However the residues of these poles contain IR 
divergent factors, meaning that via the decomposition Eq. ()3.I8j) . we can expect the Bethe-Salpeter vertex function 
Tp S and the components to also involve IR divergent factors. This will be made explicit shortly. In anticipation 
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of the cancellations necessary so that the expressions for fps and hps are indeed finite, let us define the dimensionless 
functions 



h N k = rn? A , rS + S oa - r+^-=T<& • (4-6) 



This gives the expressions 



h PS = 2^7V C / 4^ K + ] C- (4-8) 

Recall that in the chiral limit (to — > 0), Eq. (|3.26p holds. Comparing with Eq. (|4.8|) . we see that in terms of /i^ this 
is equivalent to 

-i^. (4.9) 

Given that /pg is a physical quantity and that the quark mass function, M, is IR finite, this limit tells us that 
indeed, h k should also be IR finite. Further assuming a massless pion in the chiral limit and using the definitions 
Eqs. (|4.4I4.6|) . h k in the chiral limit can be written as 

h « m ±j° —J—qn = -_Lr&, (4.io) 

which would mean that the normalized Bcthc-Salpeter vertex function component Tff obeys 

rfir^t-r-. (4-n) 

JPS 

This result agrees explicitly with the standard result (see, for example, Ref. [35j). Additionally, it confirms that the 
Bethe-Salpeter vertex function and its components are IR divergent in exactly the same way as the quark self-energy. 

Let us now turn to the homogeneous Bethe-Salpeter equation in the pseudoscalar channel, Eq. (|3.10p . Inserting 
the decompositions, Eq. (|2.7[) for the quark propagators and Eq. (|4.1[) for the pseudoscalar Bethe-Salpeter vertices, 
projecting out the Dirac components and performing the energy integrals one obtains four coupled homogeneous 
integral equations for the scalar functions E; (un- normalized at this stage). However, careful inspection reveals that 
remarkably, all four equations have integrands involving precisely the combinations of factors used in the definitions 
of the functions / and h, Eqs. (|4.4I4.6[) . The equations are: 



T o P = ~z9 C F j— [u 

2 J < w fe 

1 2 f dkP(p~k) [M+ + M-] 
1 ip - 7.9 G W -p^ r-^— =y-Jk, 



1 2 f dkP(p- k) fk [M+u k - Mj7 U +] 

T.9 J— — =r 



1 , f dkP(p-k) p-k [M++M-] 

The structure of the above set of equations shows that, in general, the Bethe-Salpeter vertex with its four components 
is a convolution integral of the IR divergent interaction, the quark mass functions (and associated quasiparticle 
energies) and the two putatively IR finite functions / and h. Notice that for equal quarks, T 2 vanishes. To arrive at 
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a closed set of equations for / and h, we must invert their definitions in terms of the IV Using Eqs. (|4.4I4.6[) . one 
may write 



P 2 
o 



+ OJp 



fv 



i 2 JP 



Op 



[M+w- + M~w+] 



-fv 



lp- 



[w+ + cj p ] [M p + + Mp 
[M+w- + M-w+] 



r 



2p" 



p 2 [M+ - M-] 



- w„ 



(4.13) 



Further, using the integral expression for the dressing function A, Eq. (|2.11[) . one can rewrite the expression for 
G) = Auj such that 



[ W P 



ctk F(p- k) p-k [uJpOJ k + uj p uj k 



[ujp + i 



(4.14) 



Substituting this integral expression and the integral forms for the Tj, Eq. (I4.12j) . into Eq. (I4.13|) one finds the final 
coupled equations for / and h: 



hp = 



fv 



P 2 



1 



tip 



2 



9 2 C F 



dkF(p- k) 



u k u k 



[M+u- + M-u+] 



V 1 M p 

dkF(p- k) 



M„ 



Ml 



[M +w- + M-w+] 



P-k [Mjuj, - M-U+] [M+ - M p 



p-k [u£+u k ] 
[MZuZ+MZw£ 



-fv 



■k [w+c 



"P -,2 

p 



P-k [w+w fc +uj p w+] 



(4.15) 



In both equations, the divergence (arising either from the 5-function charge constraint proportional to the divergent 
constant C or from the infrared singular 1 /q 4 term in F) occurs when p — k. It is straightforward to show that in 
both cases, the remaining factors of the integrand cancel. Thus, the above integral equations are explicitly infrared 
finite and so are the functions / and h, as demanded by the definitions of fps and hps- Notice that if the meson were 
anything but a color- singlet, the cancellation of the divergences would not occur, just as in the case for the heavy 
quark limit [231—2 5j . Further, we see that the only occurrence of Pq = Mp s lies in the coupling of / to the /i-equation. 
That the system of four IR divergent components of the general pseudoscalar Bethe-Salpeter vertex decomposition 
reduces to a system of two finite functions is analogous to the case of the quark propagator, where the two IR divergent 
dressing functions (A and B) conspire to form the finite gap equation for the mass function, M. 

In the special case where the quarks have equal mass, the above equations reduce significantly, although their 
general form remains the same. The equal mass equations read 



f'p hp 



dkF(p- k) 



dkF(p- k) 



LU k 



h k - h 



p-k 
P 



P 


k^ 


- M p M k 




k 2 








(4.16) 



Comparing the first equation with the gap equation for the quark mass function, Eq. (|2.13l) , one sees that the integral 
has the same form but with M replaced with h. Thus, it is clear that in the chiral limit, h is proportional to M (as 
shown in the discussion leading up to Eq. (14.91) ) and Pq = Mp s = 0, showing that the pion is massless in the chiral 
limit, as it should be. 

It is worth pointing out a major difference between the Bethe-Salpeter equation above and studies performed within 
the covariant approach, e.g., Refs. [37l - l39j . In the covariant gauge approach, one works in Euclidean space and since 
the bound state total momentum is timelike, it is necessary to make an extension of the quark propagator to complex 
Euclidean momenta (usually numerically). This extension into the complex plane has the consequence that possible 
singularities must be treated carefully or avoided altogether (although numerical techniques to deal with such issues 
have been developed, e.g., Refs. [40h43| ) . For mesons of roughly equal quark masses this is not a major issue, however, 
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for asymmetric systems such as heavy-light mesons this becomes challenging. In the Coulomb gauge rest frame (and 
Minkowski space) equations presented above, only the quark mass function (and quasiparticle energy) enters and is 
always evaluated at spatial momenta, rendering the issue of continuation to timelike momenta irrelevant (effectively, 
the energy integrals performed analytically in the derivation of the equations deal with this issue from the outset). It 
will be seen that heavy-light systems are no more difficult to treat than equal quark systems. 

In order to calculate fps and hps, we must normalize the functions / and h. As a first derivation, let us consider 
the equal mass case where 

fps = 2*N C [ (4.17) 
J w fc 

However, in the chiral limit, Eq. (|4.9p holds such that we can write 

M k n ^2if PS h%, (4.18) 
which would mean that the normalized functions must obey 

i ^ - ANc f **. f » h » (4.i9) 

in this limit, such that the GMOR relation, Eq. (|3.25[) . is satisfied. We see that in this restricted limit, the normal- 
ization can be written in terms of our finite functions / and h. More properly, the canonical normalization of the 
Bethe-Salpeter vertex function reads (in the Minkowski space rest frame and with trace over color, flavor and Dirac 
matrices) [35[ 



*8*P° = Tr f , c , d Jdk{rf s (k:~P°)^^ 

(4.20) 

Evaluating the flavor (recall that we have Nf = 2 and consider flavor nonsinglet mesons) and color traces, and using 



one obtains 

1 = -^7 Tr <* / dkf N ps (k;~P°)W+(k+) {j Q W+(k+)T^ s (k;P ) - I%(£; P*)W- q (k-)f) W^k~). (4.22) 

The (Dirac structure of the) conjugate Bethe-Salpeter vertex function is explicitly given in the rest frame by 

T N (k; -P°) = CT NT (-k; -P )^ 1 = 7 5 [r& P° 7 °rf & - - P^j-kn] . (4.23) 

After evaluating the Dirac trace and performing the energy integrals, one finds after some effort that the normalization 
condition can be written in terms of the finite functions / and h: 

1 = -4AT / d%f " H " t M * + + M *1 (4 24) 

which reduces explicitly for equal mass quarks to the earlier form Eq. (|4.19j) . Given an arbitrarily normalized solution 
to the Bethe-Salpeter equation, Eq. (|4.15l) . (we shall use the condition h(p — > 0) = 1 for the numerical analysis later 
on) it is useful to write 

This then results in the normalization condition 

N = AN C / — 4^ r T U ~TT (4-26 
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and the redefinition of fps and hps to 



_2N C f dk [M+ + M-] 
Jps - -rr \ + _ r + , rWfc, (4.27) 



IPS 



2N C f dk 
~JT J oj+ojI 



(4.28) 



Recall that it is known that the Coulomb gauge quark propagator under this truncation reduces to its heavy quark 
limit (in the absence of pure Yang-Mills theory corrections) [30(. The difference between the mass function and the 
bare mass becomes smaller in the heavy quark limit, i.e., Mk - m/, -> as m/, -> oo (we will use m/j to denote the 
heavy quark mass). For a heavy-light system, with m + = mj, — > oo and assuming that the above integrals for J\f, fps 
and hps converge for large k (i.e., that the convergence is driven by the large fc behavior of the functions fk and hk) 
then we can make the replacement Lu k — MjT = 3> uJ k , MZ within the integrals. This then gives 

, 2 m+=m h 4iV c f dkf k h k m+=m h 2N C [3k m+=m h 2N c ftik 

JV — > / _ r _ , : ,_-| , Jps — > 77 / -^r/fci h PS — > -77- / -^h k . (4.29) 

m h J u k [co k +M k \ Mm h J uj k At J cu k 

Performing the same substitution for Eq. (|4.15[) gives 

hp m+ ^- + 0(l/m h ), f p m+ -^ h h p + 0(l/m h ), (4.30) 



in 



h 



which tells us that (asymptotically, as mh — > oo) Pq = Mp s = m\ as we would expect. We notice though, that the 
leading order heavy quark limit merely tells us that f p = h pi but not what f. p or h p actually are. With the above 
limits for Mps, fps and hps, we see that the relationship Eq. (|3.24j) is explicitly fulfilled. The most important result 
of this analysis though, is that as mh — > oo (assuming the integrals are well-defined) and given M ~ l/y/mZ, 



fpsVMps - const. (4.31) 

This result explicitly agrees with the analysis of heavy quark effective theory at leading order [44j . 

At this stage, it is appropriate to make a comparison of the formalism presented here with similar works on the 
Coulomb gauge Bethe-Salpeter equation. As mentioned earlier, the chiral limit pseudoscalar meson (i.e., the pion) 
was investigated in Ref. [26j. An extension to the case of equal, but finite mass quarks was performed in, e.g., 
Refs. [H,|4^, [4(|. The equations therein have subtle differences to those presented here, over and above the obvious 
extension to arbitrary quark masses. Since the referenced works dealt with equal mass quarks, the vertex component 
T2 was absent from the outset although they explicitly found that the system can be decomposed in terms of two 
finite functions. The coupled equations for the two finite functions (/ and h) presented here are, from the point of 
view of canceling the IR divergences, somewhat more explicit than the earlier versions since the IR divergent factors 
(cD in our notation) inherent to the quark self-energy have been completely eliminated. This elimination also extends 
to the equations for fps, hps and the normalization. The success of making this cancellation explicit arises from 
utilizing a different definition of the two functions, Eq. (|4.6p . designed with the arbitrary quark mass expressions for 
fps and hps in mind. This is in distinction to the original decomposition of Refs. [26l. |28|. whose motivation was 
to make the Bcthc-Salpcter equation tractable. Despite the differences in the decomposition, both sets of functions 
yield finite quantities. 

V. VECTOR MESONS UNDER TRUNCATION 

Let us now consider vector mesons. Although the expressions are considerably more involved than those discussed 
in the previous section, they follow the same pattern. The homogeneous Bethe-Salpeter equation retains its form as 
for the pseudoscalar case, Eq. (|3.10l) . 

r a v MP Q ) °= v -*g 2 c F J ^fcF(p-fc) 7 ^+(fc + )r^(fc;P°)w^(fc-) 7 , (5.1) 

but the general decomposition of the normalized vector Bethe-Salpeter vertex function (in the rest frame and omitting 
the flavor dependence) now reads 

rm p°) = {Y K P + p° 7 °rf p + t?t» + p°7°rp-rf ] 

+p l + r& + p° 7 (it P - rg) + rpvl + p^tpT"] } , (5.2) 
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where it is recognized that the vector meson at resonance is transverse to the total four-momentum P and so, can only 
have spatial components. The vector meson leptonic decay constant can be defined as (trace over Dirac matrices) [38[ 



fv = - 



^-Tr d J dkfW+(k + )T^(k;P )W^ q (k-) 



(5.3) 



and the normalization defined as (trace over color, flavor and Dirac matrices) 



6i5 ab P° = Tr />C;d / dkl Yyl (k; -P u ) 



Y™ (k; P°)W-(k-) + C (k; -P'WUk^T™ (k- P°) 



dP° 



dP° 



In both the latter equations, a factor of three arises from the (transverse) polarization states of the vector meson. 
With the definition of the conjugate amplitude as before: 



T Vl {k--P») = CT»t{-k--P»)C 



(5.5) 



one has 



r^(fc; -p°) =y -r& + p° 7 °rft + ^ k + p^-kr 



:>>k 



(5.6) 



Starting with the leptonic decay constant, fy, defined in Eq. (15.3[) . after inserting the decompositions for the 
Bethe-Salpeter vertex, Eq. (|5.2p . and the quark propagator, Eq. (|2.7[) . evaluating the trace over Dirac matrices and 
performing the energy integrals, one has the expression 



fv = ~ 



2iN r 



dk 



3MvJ W ^ fc -[P 2 -K + +^) 2 



[K + M k] 



(^ok [K"k + M k -u£] + rf fe fc 2 [u,+ - ^]) 



^Ok + k 2 ^Qk 



M"k + M k -"t) rf fe fc 2 [w+ + W k 



2rf fc [M+u> k + M k u+] + 2lf fc fc 2 [cj+ + + 



1 lk ft 1 7fe 



[K^k + M k ut] + lf fe fc 2 [ w + - W fc ]] } . 

(5.7) 



Interestingly, we note the appearance of the combinations of functions Qq and , Eq. (|4.4[) , defined for the pseu- 
doscalar meson. Introducing the further two combinations 



r 



'2 P JV 
6fe 



fc^r 



1 lfe K 1 7fe 



[M^k + M fc -o;+] - lf fe P [c+ + a;"] , 

[M+ Wfc - + M fc " W +] + if fc P [ w + - W j] , 



(5.. 



we can write 



fv = 



2iN r 



dk 



3M V J u+uJ k [p2-(Q++Q k y< 

M +uj 



X \ (Wfc + ^fc ) 



[M+ + M fe -] 
"Vofe r + ^i - V2fe 



K" + w fe ] 



- P 2 [2Q& + Q&] 



(5.9) 



We identify the appearance of h k (defined in Eq. (|4.6[) for the pseudoscalar meson) in this equation and introduce, 
in analogy to before, the dimensionless functions 



i 



s * + ) M + u k \ nN 1 
— V3fe r ' 



m+M- 



,-.N 
Jk 



P 2 



[P 2 -(^ 



2 ] \H+"k] [ M k + M k 



(5.10) 
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where the function gjf will shortly be needed and the aim will be to see if these functions are finite. The vector meson 
leptonic decay constant can thus be written as 

fv = S / !§F k (2 ^ + w *1 ^ + t Mfe+ + M ^ 3 ^ • (5 - u) 

The above expression has an almost identical form as for the pseudoscalar case, Eq. (|4. T[) . 

Turning now to the vector meson Bcthe-Salpeter equation, Eq. (|5.1[) . after inserting the decompositions for the 
Bethe-Salpeter vertex, Eq. (|5.2j) . and quark propagator, Eq. (|2.7I) . projecting out the Dirac structure and performing 
the energy integrals, one can identify the appearance of precisely the combinations of (un-normalized) functions that 
define f k , h k , g k and j k , Eqs. (|4.6l5.10p . (the former two defined for the pseudoscalar meson) in exactly the same way 
as previously. That there are only four combinations means that there are only four independent dynamical functions 
that characterize the vector meson under this truncation. The equations for the vector meson I\ read: 

To P =\ 9 2 C F [ ^ F }P -V K+ + h k - | [M+ + M-] Jk 




1 m f <tkF(p-k) j ^ A\ [M+ + M-J j, , A 



r 2P =\<?c F [ ^F(pJ) p^_ { _ hk] 

1 2 f dkF(p-k) p-k [M+ + M fc ~] 
lsp-^5 Oj? / -qr-- 35- r „,+ ,,_ : ,,_ ,+i/fc, 



2 J *>t u k P 2 [ M k"k+ M k"t] 
1 2 /• dfcF(p - fc) p-fc [<4 + uj\ [M+ + M fe ~] , 

1 2 f dkF(p-k) p-k K-<j 
J-Bp— r0 Of / q— : — r + _ +1 fffc, 

r 6p =^9 2 Cf / ^— ^— 1 - — { fof + u; fe /i fe - M+ + M fe ] j fe } 

lj, /" dkF(p-k) 1 / 3A\/[M++M fc - 



P 




r 7p =-<? / ,.+ .- — ( 1 - — ) \ r..+ , ..-i /* - ^ > . ( 5 - 12 ) 



where we write 



p-k 2 
p 2 k 2 



(5.13) 



Assuming that the functions Jk, h kl gk and jk are finite (as we shall shortly demonstrate), the vector meson Bethe- 
Salpeter vertex IV (and its components, I\) is given as a convolution integral involving the IR divergent interaction 
and the quark mass functions, just as for the pseudoscalar meson. Further, we notice that if we artificially set A = 0, 
the first four equations are identical to those in the pseudoscalar case. In the special case of equal mass quarks, r 2 
and r 5 vanish. Inverting the definitions of g k and j k , Eq. (I5.10[) . one has 

Pq9p - [&% + &p] [w+ + u~] j p =(r 0p + p 2 T 6p ) [M+uj- + MpUj+] - T 4p p 2 [uj+ + u)~] , 
•T.+ j_,7i-1 lM+,,,- 4- m-,.,+ 



r +^ -i 9p-Jp =( r ip-P r 7 P ) r . f+ - Af _i +r 5p p f— q— 1 (5.14) 

[up + uj p \ [Mp + M p \ [Mp + Mp \ 

in addition to the analogous expressions defined as for the pseudoscalar case, Eq. (|4.13[) . Collecting together the 



15 



definitions, along with Eq. (|4.14[) for the factors involving ui, this results in the following four equations: 
h„ = 



P 2 




1 2 /• dkF(p-k) 
Jp + 7:9 C F q-^ 

l J uTu,. 



A 



[M+u- + M-u+] 



[u+ + up] [Mp 

f k [M+w^ - M^u+] [M+ - M-] p- k [w+Wfc + u~u+] \ 




p -2 r + i -i 
p [up + uj p \ 



fp = h P + o3 C 



1 A\ [M+ + M fc -] 

2 / k + + ^1 



[M+u- + M-cj+] 



[<4" 



- U v 



p-k [u+ 



1 



Jp = 



P 2 
o 




2 5p + 2# ^ 



p 2 [wp + w p ] 
dk F(p — k) 



J k ^k 



[Mt + M-] ]k + A[u++u- k ] h k 



[M+u- + M-u+ 



+jk 



P-kH+^k] [Mk+M^ 



[u++Up][M+u^ + M,;up h ■-- 



fk [u+u k +u p u+ 



9 P =j P + 



dkF(p- k) 



p-k 2 
p 2 k 2 



P 

9k + A 



[M+w- + M-w+] 



[M+ 



P'fc Wi 



1 K 



[M+ Wi 7 + M" W +] [M+ + Mp 



9p 



P-k [u+u k + u p u+ 



P 



[up 



(5.15) 



As before for the pseudoscalar mesons, the singularities arising from F when p 
functions are indeed finite. In the equal mass case, the equations reduce to 



k cancel, showing that all four 



~J P + 




p-k t AM k 

2 u k 



"P ->2 
P 



Jk 



f P = h P + 2 9 Cf 



J P =^9p + \9 2 C F 



X -2 



A\ M p M k 



u k | \ * / ^ k 

dkF(p~k) j p-k 2 M p M k 
p 2 k 2 u p u k 



Uk 



Jk T o J* Jp ^2 
^k P 

p-k .p-k 

Jk ~ Jp~/2 
UpU k p 



-Jk 



1 2 f dkF(p-k) p-k 2 Up 
9p =J P + o.9 C F / <^ ^2T7 — 

2 J u k \pk 2 Uk 



9k 



p-k 
9p~a 
P 



A *pMk fk 



ur 



A M„ 

"o" 9k 

2 w fe 



A^. fc 

Up 



(5.16) 



The above sets of coupled equations have a rather illuminating structure. One sees that outside the integrals, the 
functions (/, h) and (g, j) appear pairwise, each with only one occurrence of Pg = My. Inside the integral, the pairing 
is (/, g) and (h,j). This pairwise connection will be visible in the numerical solutions. Within the integrals, the pairing 
comes with an associated factor of A that serves to couple the two sets of equations for the pairings (/, h) and (g,j) 
(hence the decoupling when A is artificially set to zero). A, as defined in Eq. (|5.13[) , is dependent only on the angle 
between the spatial momenta p and fc. That when an angular dependence is suppressed, the equations for the vector 
meson reduce to those of the pseudoscalar has an obvious interpretation: the emergence of the quantum mechanical 
description of the vector meson as being a total angular momentum excitation of the groundstate (pseudoscalar) 
meson, even in the context of the relativistic (though noncovariant) Bethe-Salpeter framework used here. 

The normalization of the vector meson Bethe-Salpeter vertex function can now be discussed. Taking Eq. (|5.4p . 
inserting the vector meson vertex decomposition, Eq. (|5.2p . its conjugate, Eq. (|5.6[) , along with Eq. (|2.7p for the quark 
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propagator and treating as for the pseudoscalar case, the result is 

3 J lo+lo- L J [M+u k + M fe "^] V 7 

or, for equal mass quarks 

l = -^/4[2/W + ff fif], (5.18) 
" fc 

and where again, only the finite combinations fh, hk, gk and jk defined in Eqs. (|4.6I5.10[) appear. Just as for the 
leptonic decay constants and the Bethe-Salpeter equations, the vector meson normalization has the same structure as 
for the pseudoscalar case, Eq. (|4.24|) . 

Finally, let us discuss the heavy quark limit for the vector meson. This proceeds as for the pseudoscalar case. With 
the normalization convention that 

(h p and f p as in Eq. (|4.25[) and with — > 0) = 1), we have from the normalization that (m + = rrih — > oo as before) 

J* — > \ / - r - , ^-1 2 /* ft * + 9k3k] , (5.20) 

or Af ~ 1/y/mh- The vector meson leptonic decay constant, Eq. (|5.1ip . becomes 

* -> 3A^y^ (2/ifc+jfc) - (5 - 2i) 

The vector meson Bethe-Salpeter equation, Eq. ()5. 15[) . gives us the further information that f p — h p , j p = g p and 
My = m/j in the heavy quark limit, such that as — > oo, 

fv\/Mv~ const. (5.22) 

Again, this agrees with the known heavy quark limit at leading order fiij . There is one further result that we should 
check and this concerns the ratio fy / fps ■ The heavy quark effective theory result at leading order is 44] 

p_ = 1 _ ^pl (5 . 23) 

J PS 37T 

(a s being evaluated at the scale m/j). In our case, we do not include a s contributions, so we should see that for 
heavy-light systems as m + = — > oo 



P- m+ -^> h 1. (5.24) 
J PS 



VI. NUMERICAL IMPLEMENTATION 

Before presenting the results, it is worthwhile discussing the numerical implementation of the equations. While 
the infrared singularities of the interaction do cancel formally, they present a formidable challenge for the numerical 
analysis. The techniques used are based on those employed to solve the gap equation for the quark mass function 
M k , Eq. (|2~T3) . explained in detail in Ref. 

As a concrete example, let us consider the equal quark mass pseudoscalar meson Bethe-Salpeter equation for h p , 
Eq. (|4.16p . which reads 



p ° f ,1 i r [ dkF(p-k) f fk\ 

^f P + -9 C F J \h k -h p y\ (6.1) 
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(all the other Bethe-Salpeter equations follow the same pattern). In conjunction with the partner equation for f p , 
the above is one of a homogeneous set of equations, valid only at a particular value of Pq = Mp s . The standard 
technique would be to first perform the angular integral involving F(p — k) (in this case, the integral can be performed 
analytically). Then, the radial integral would be set up on some discrete grid of k 2 and p 2 values which results in 
a matrix form for the coupled equations. One can then use standard matrix methods to determine the eigenvalues 
(the groundstate meson mass corresponds to the value of Pq where the highest eigenvalue is equal to one) and 
eigenfunctions. However, in the above case the matrix would involve cancellations of singularities (P is divergent 
when k = p) between the various terms. To circumvent this would in principle require derivative information about 
the functions or an infrared regularization of F such as a fictitious mass term. Further, when k is in the vicinity of 
p, the matrix elements become significantly enhanced. Keeping control over the numerical fidelity in such a situation 
would become extremely difficult. 

As pointed out earlier, the above equation in the chiral limit is very similar to that of the quark gap equation, 
Eq. (|2.13|) . whose numerical evaluation is well understood [30(. The first step is to change integration variables such 
that the radial momentum flows through F . Denoting q = p — k, this then gives 

h Po , 1 2 r f dkF k f p-q\ 

hp = -r^fp + t;9 C F / < h q - hp^ \ . (6.2) 



Because of the infrared singularity, it is useful to rewrite this equation in the form 

X p + Ip 

TT7T 



K = -TTwT. (6-3) 



P 2 , 1 „ f dkF k , 1 , f dkF k p-q , s 

X v = Ip = ~ 9 2 C F / — ^ V I 2 P = -g 2 C F / — (6.4) 



where 



^ 2 p V r * J UJ q " T J UJ q p- 

This manipulation substitutes the difference of two infrared singular integrals for their ratio, which is numerically less 
prone to instabilities. Along with the corresponding equation for f p (manipulated in the same way), we now introduce 
a fictitious eigenvalue A(Pq ) to the left hand sides of both equations, such that we have the form (similarly for f p ) 

n X r) -(- p} 

KP 2 ,)h P = -f^jf- (6.5) 

Taking a starting guess for the functions h p and f p with the boundary condition h(p 2 — > 0) = 1, we construct the right 
hand sides of the equations. This is doing by interpolating or extrapolating h q and f q within the angular integrals 
as required, using standard techniques. As will be seen, the functions h p and f p are smooth such that this is a 
numerically safe procedure. The radial integral is performed over a discrete grid of values for k 2 , regulated using an 
infrared cutoff e and an ultraviolet cutoff A (both with dimensions of [mass] 2 ). The derived value for h p at the lowest 
p 2 is now compared with the boundary condition (unity) to give the eigenvalue A. Both derived functions h p and 
f p are then scaled by this value such that the boundary condition is fulfilled. The coupled set of equations is then 
iterated until convergence of both the functions and the eigenvalue is achieved. Essentially, this procedure is identical 
to solving an inhomogeneous equation iteratively, but with the lowest grid point value for the function (fixed to unity 
here) swapped for the eigenvalue (which for the inhomogeneous equation would be trivially absent). The value of Pq 
for which A = 1 is then found. 

It is worth pointing out that the numerical procedure described above is certainly not without its limitations. Both 
integrals I p and I 2 diverge as as e — > 0, such that the iteration steps are damped (typically by a factor of 

10 3 ) and one must use many steps to achieve a particular tolerance. However, this convergence is stable and easy 
to implement. If one were to use the difference of the two integrals rather than their ratio, errors in the unknown 
functions are amplified by such factors and the convergence is highly unstable - one would require an extremely precise 
starting guess for the functions in order for the iteration to work. 

Once the functions h p and f p (and in the case of the vector, j p and g p too) are found, it is straightforward to obtain 
the normalization and the leptonic decay constant. In the case of the pseudoscalar meson, hps can also be found in 
order to verify the relation Eq. (I3.24|) . 
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VII. RESULTS 



Aside from the quark masses, there is in principle only one free parameter in the truncation scheme studied in this 
work and this is the string tension a used in the input function F. By setting a = 1, all quantities are then (and 
unless otherwise stated) expressed in appropriate units of a. To compare with physical quantities, we also give an 
estimate of the magnitude using the range \fo ~ y/aw — y/^ a w ~ 440 — 762MeV as discussed in Sec. II. We use 
the following numerical parameters throughout this section: the infrared cutoff, e = 10~ 6 , and the ultraviolet cutoff, 
A = 10 6 . Concerning the value of e: in Ref. [3(| it is seen that the quark mass function converges for e < 10~ 6 . 
We have checked the results here against e = 10~ 5 (given the infrared singularities involved, lower values of e result 
in significantly more computational effort), and the difference is typically < 10~ 3 (< 3 x 1CP 3 ) in the pseudoscalar 
(vector) meson masses and <^ 1CP 4 (< 3 x 10 -3 ) in the pseudoscalar (vector) meson leptonic decay constants. With 
the typical range of a values, this tolerance would correspond to ~ lMeV. However, notice that for the quark mass 
function, the difference between results is much smaller for e — 10~ 6 — 10~ 7 than between e — 1CP 5 — 1CP 6 [30| . 
Thus, we consider the results presented here to be convergent within reasonable limits. The value of A is chosen such 
that the heavy quark masses (up to m = 10 3 ) can be accommodated: with the interaction Eq. (|2.3|) . all results are 
UV-finitc. 

Let us begin with the case of (equal mass) chiral quarks. Solving Eq. (|2.13[) for the quark mass function (as in 
Ref. do[) and inserting into Eq. ([QUI) , the chiral condensate is found to be -0.0123, which with ■ s /a = 440 - 762MeV 
gives <qq> = —(102 — 176MeV) 3 . While the quark condensate is not a physical quantity, this value is considerably 
lower than the commonly accepted <qq>^ — (250MeV) 3 . The result here is consistent with previous similar studies, 
e.g., Refs. [13, [29| . Turning to the pseudoscalar meson, this is a special case in that we know analytically that Mps = 0. 
The pseudoscalar leptonic decay constant is found to be fps = 0.0260, which corresponds to fps ~ 11 — 20MeV with 



27]. Comparing with the physical 
36j, we see that the pseudoscalar 



the above range for </cr. Again this result agrees with previous studies, e.g., Ref. 
result for the pion (with our conventions for the normalization), /„. = 92.2MeV 
leptonic decay constant for light quarks is almost an order of magnitude too small. Clearly, while chiral symmetry 
is dynamically broken by the infrared enhanced interaction, the absence of the spatial components of the interaction 
makes a big difference to the amount of symmetry breaking. The inclusion of more detailed interactions has been 
recently studied within this context in Ref. (29j . 

Turning to the chiral quark vector meson, we find that My = 1.553, or My — 683— 1184MeV. This has the correct 
magnitude when compared to the physical result, M p — 775.5MeV |4 7| . The leptonic decay constant is found to be 
fv = 0.250, or fy — 110 — 190MeV, as compared to f p = 153MeV [38j (with our normalization convention). The 
calculated vector meson leptonic decay constant is thus also comparable with the magnitude of the physical result, 
unlike its pseudoscalar counterpart. 

Staying with the chiral mesons, the normalized vertex functions h N , f N and for the vector, g N and j N , are plotted 
in Fig. |31 Recall that these functions are constructed so as to be free from singularities and this is explicitly verified 
(and has been checked for all combinations of quark masses). The functions are smooth and well-behaved such that 
the numerical implementation of the angular integrals is indeed justified (see the previous section). In the chiral 
limit, where Mps — 0, the pseudoscalar function h should be identical to the chiral quark mass function, M, up 
to a constant factor and this has been verified (and serves as a useful numerical check). The vector functions h N 
and j N have a similar form, but with minor differences. The functions / (pseudoscalar and vector) and g N exhibit 
a bump at around x — k 2 ~ 0.2 — 0.5, although this does not seem to have any significance. In the vector case, 
the pairwise behavior of the functions (/, g) and (h, j) is clearly visible in the IR. As mentioned in the discussion 
following Eq. (|5.16p . these pairings occur in the integral expressions for the vector meson Bcthe-Salpeter equation: 
the pairs have the same infrared limit, but at finite momenta j N and g N are generally smaller and shifted slightly 
towards infrared momenta. 

Let us now consider mesons with equal mass quarks. The meson masses are plotted in Fig. [4] as a function of the 
quark mass, m. One can see that the pseudoscalar meson mass goes approximately like y/m in the chiral limit, a clear 
signal for the dynamical breaking of chiral symmetry and the Goldstone boson nature of the pion. For large quark 
masses, both meson masses approximate to 2m and are degenerate. We have verified that Eq. (|3.24[) is numerically 
satisfied (this again serves as a useful numerical check). These results are in explicit agreement with those presented 
in Ref. [4^]. The binding energy, defined as M PS — 2m (similarly for the vector meson), is also plotted. This provides 
a more detailed picture of the interaction, since the meson mass is dominated by large quark masses in that limit. 
One can see that the binding energy is not trivial and decreases for larger quark masses. The degeneracy of the meson 
masses for heavy quarks is also clearly visible from the binding energy. That the pseudoscalar and vector meson masses 
are almost equal for equal quark masses above 0(1), means that the mass splitting of charmonium and bottomonium 
is not present. The leptonic decay constants for the equal quark mass mesons are plotted in Fig. [5] Unlike the meson 
masses, the leptonic decay constants are only degenerate in the heavy quark limit. For bottomonium, recent lattice 
results give fps — 472MeV 48] (with our conventions); comparing with our results, we get (for the heaviest quark 
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FIG. 3: [left panel] Pseudoscalar and [right panel] vector normalized vertex functions with (equal) chiral quarks, plotted as a 
function of x = k 2 . All dimensionfull quantities are in appropriate units of the string tension, a. See text for details. 
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FIG. 4: Pseudoscalar and vector meson masses [left panel] and binding energies [right panel] with equal mass quarks, plotted 
as a function of the quark mass. All dimensionfull quantities are in appropriate units of the string tension, a. See text for 
details. 



mass) at most fps ~ 0.48, or fpg ~ 211 — 366MeV. The vector meson result is similar, meaning that both leptonic 
decay constants are significantly too small, even taking into account the range of values for a. In the equal quark mass 
case for large quark masses, the pseudoscalar and vector normalizations become degenerate and go to zero; further, 
all vertex functions become degenerate. 

Finally, let us consider the case where one quark has fixed mass. The meson masses are plotted in Fig. [S] for the 
case where one quark is fixed to zero bare mass (i.e., m~ = 0, m + — to). In the limit as the other quark mass (m) 
vanishes, Mps ~ \/m (signaling dynamical chiral symmetry breaking), whereas for large masses, Mp$ ~ My ~ to 
as one would expect. Again, Eq. (|3.24|) has been numerically verified. As previously, defining a binding energy as 
Mps — m (similarly for the vector) reveals the interaction in slightly more detail and this is also shown in Fig. [5] 
The meson masses are degenerate above to ~ O(10) and the mass-splitting between states is again too small in 
this limit. Unlike the case of equal mass quarks we see that the binding energy does not decrease for large quark 
masses. Turning to the leptonic decay constants, plotted in Fig. [7J the results are again consistently too small. For 
example, fp = 144MeV [4^, [H3l (with our convention). Using the input value \fa = ^fow and a quark mass of 10, 
corresponding to a bottom quark of to = AAGeV in these units, one has fps = 0.168 ~ 74MeV, which is clearly too 
low. On the other hand, taking ^/a = y/3aw and a quark mass of 4.64, corresponding to to = 3.5GeV, one has the 
result fps = 0.190 ~ 145McV: at first glance this looks like a nice result, but given that the bottom quark mass (in 
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FIG. 5: Pseudoscalar and vector meson leptonic decay constants with equal mass quarks, plotted as a function of the quark 
mass. All dimensionfull quantities are in appropriate units of the string tension, a. See text for details. 
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FIG. 6: Pseudoscalar and vector meson masses [left panel] and binding energies [right panel] with one fixed chiral quark, plotted 
as a function of the other quark mass. All dimensionfull quantities are in appropriate units of the string tension, a. See text 
for details. 



physical units) is realistically too low and that the leptonic decay constant decreases from this point onwards as the 
quark mass increases (see Fig. [7]), the result would be also too small. However, the real interest lies in the heavy quark 
asymptotic behavior of the leptonic decay constants. In both cases, we should see that the leptonic decay constant 
decreases as one over the square root of the meson mass and this is explicitly verified in the right panel of Fig. [7J 
where it is seen that Eqs. (|4.31l5.22j) are fulfilled. Further, it is clear that Eq. (|5.24p is satisfied. Looking at the 
vertex functions in the most extreme case (m = 10 3 ), plotted in Fig. [5J one sees that the functions are all identical. 
The heavy quark analysis told us that f N = h and j = g ; seemingly the pairing of h with j N and f N with 
g for the vector meson, seen for the equal chiral quark case, also happens here. Thus, the numerical solution to the 
truncated system does indeed reflect the heavy quark symmetry. 
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FIG. 7: Pseudoscalar and vector meson leptonic decay constants [left panel] and fpsVMps, fvVMy [right panel] with one 
fixed chiral quark, plotted as a function of the other quark mass. All dimensionfull quantities are in appropriate units of the 
string tension, a. See text for details. 
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FIG. 8: [left panel] Pseudoscalar and [right panel] vector normalized vertex functions with one chiral quark and one quark of 
mass 10 3 , plotted as a function of x = k 2 . All dimensionfull quantities are in appropriate units of the string tension, a. See 
text for details. 



VIII. SUMMARY AND CONCLUSIONS 



In this study, a leading order truncation scheme developed in the context of Coulomb gauge Dyson- Schwinger 
equations [30j has been extended to the Bethe-Salpeter equation framework for quark-antiquark pseudoscalar and 
vector meson masses and leptonic decay constants, with quark content of arbitrary mass. The leading order Coulomb 
gauge truncation scheme centers around the inclusion of an instantaneous and IR singular interaction corresponding to 
a pure linear rising potential, supplemented by a term arising from the consideration of total color charge conservation. 
As was briefly discussed, this results in a particular form for the quark gap equation (the Adler-Davis gap equation 
f27l ]) for the static quark mass function and for which the divergences cancel and dynamical chiral symmetry breaking 
is observed. The quark proper two-point function and the self-energy are divergent, reflecting the physical picture 
that one would require an infinite amount of energy to create a colored object from the vacuum. Further, under 
this truncation scheme, a connection to the Coulomb gauge rest frame heavy quark limit (in the absence of pure 
Yang-Mills corrections) |23[ can be established. 

The connection of the quark propagator to the quark-antiquark (flavor nonsinglet) Bethe-Salpeter eq uation is 
well understood via the axialvector Ward-Takahashi identity (AXWTI) . Following the approach of Ref. [34j (written 
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within the context of the covariant formalism), the necessary framework was briefly reviewed and adapted to the 
(noncovariant) rest frame Coulomb gauge formalism studied here. The importance and usefulness of the AXWTI 
is primarily in that it guarantees that the pion emerges as the massless Goldstone of dynamical chiral symmetry 
breaking in the chiral limit. Indeed, earlier studies in Coulomb gauge were focused on this aspect [26|, [28[. However, 
the AXWTI also applies to arbitrary quark masses, including heavy quarks and allowing a simultaneous analysis of 
both the chiral and heavy quark limits. 

The necessary formalism for calculating the pseudoscalar and vector meson masses and leptonic decay constants 
was derived. Given the IR singular interaction along with the associated divergence arising from the consideration of 
total color charge conservation, the overriding concern was the elimination of the IR divergences within the equations. 
It was seen that this can be explicitly achieved for arbitrary quark mass content and that the pseudoscalar meson 
is described, within this Coulomb gauge truncation scheme, by two expressly finite functions; the vector meson 
correspondingly described by four finite functions. In contrast, the Bethe-Salpeter vertices themselves are manifestly 
divergent. This situation is directly analogous to the case of the quark with a finite mass function though divergent 
self-energy. Such cancellations have also been observed in the Coulomb gauge rest frame heavy quark limit for mesons 
[23j|, baryons [24l ]. and in the four-point quark Green's functions [25I ]. The connection to the heavy quark limit was 
presented and it was shown that the leading order results of heavy quark effective theory [44| are reproduced. 

The Bethe-Salpeter equation reduces to a particularly intriguing form within this truncation scheme. As discussed, 
the expressions, Eqs. (|4. 1515.151) . are a generalization of earlier Coulomb gauge studies of the pseudoscalar meson 
[26l [28j to arbitrary quark masses and extended to include vector mesons. Given the fact that the Bethe-Salpeter 
equation is constructed from the quark proper two-point function via the AXWTI, it is no surprise that the equations 
for the pseudoscalar meson have a very similar structure to the quark gap equation (27l . [30j , and directly exhibits the 
chiral limit properties that one would expect. The dependence on the bound state energy arises in simple fashion and 
the vector meson appears as an angular momentum excitation of the pseudoscalar. In contrast to covariant gauge 
studies of the Bethe-Salpeter equation (e.g., Refs. 37-39]) there are no complications, at least within this leading 
order truncation scheme, arising from the necessity to evaluate the quark propagator dressing functions at complex 
Euclidean momenta (i.e., the timelike regime) and dealing numerically with possible non-analytic structures (such as 
in Refs. [40l - [43j ] ^ . This allows for an analysis, both analytic and numerical, of mesons with arbitrarily asymmetric 
quark content. 

The derived Bethe-Salpeter equations were solved numerically and results for the pseudoscalar and vector meson 
masses and their respective leptonic decay constants were presented for arbitrary quark masses. It was explicitly 
verified that the IR singularities cancel. Clearly visible were both the qualitative pattern of dynamical chiral symmetry 
breaking and the leading order heavy quark limit. In this respect, the leading order truncation, despite its lack of 
sophistication is rather successful. Quantitatively, the degree of dynamical mass generation, the mass splitting between 
the pseudoscalar and vector meson masses, and the leptonic decay constants are all too small (although the vector 
meson mass and the leptonic decay constant in the chiral limit do have roughly the correct order of magnitude). 
Having focused on the cancellation of infrared divergences in this study, the coupling of the quarks to the spatial 
components of the gluon field and the gluon self-interaction have been neglected within the context of the leading 
order truncation scheme. Clearly, a quantitative comparison to the physical spectrum should include such interactions 
and this we propose to do in future work. 
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